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Abstract 

This paper is devoted to the study of the stability of limit cycles of a nonlinear 
delay differential equation with a distributed delay. The equation arises from a model 
of population dynamics describing the evolution of a pluripotent stem cells population. 
We study the local asymptotic stability of the unique nontrivial equilibrium of the delay 
equation and we show that its stability can be lost through a Hopf bifurcation. We 
then investigate the stability of the limit cycles yielded by the bifurcation using the 
normal form theory and the center manifold theorem. We illustrate our results with 
some numerics. 

Keywords: Delay differential equations, distributed delay, Hopf bifurcation, stability, limit 
cycles, normal form, center manifold, blood production system, hematopoietic stem cells. 
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1 Introduction 



This paper is devoted to the analysis of the nonlinear delay differential equation 



x'(t) =-(S + (3(x{t)))x(t) + - ( /3(x(t - s))x(t - s)ds. 



(1) 



This equation arises from a model of pluripotent hematopoietic stem cells dynamics, that is 
stem cells at the root of the blood production process [U [2] • It describes the fact that the cell 
density evolves according to mortality and cell division. One may stress that the cell density 
considered in equation ([1]) is in fact the density of resting cells, in opposition to the density 
of proliferating cells. 

The distinction between these two stages of the cell cycle is now widely accepted. We 
can cite, for example, the works of Burns and Tannock [3] on the existence of a resting 
phase — also called Go-phase — in the cell cycle. This phase is a quiescent stage in the cell 
development, contrary to the proliferating phase which represents the active part of the cell 
cycle: it is composed of the main stages of the cell development (e.g. DNA synthesis) and 
ends at mitosis with the cell division. Thus, proliferating pluripotent hematopoietic stem 
cells are committed to divide and give birth to two daughter cells which immediately enter 
the resting phase and complete the cycle. 

Mathematical models describing the dynamics of hematopoietic stem cells population have 
been studied since the end of the seventies and the works of Mackey [H [5] . For the reader 
interested in this topic, we mention the review articles by Haurie et al. [B] and Mackey et 
al. [TJ, and the references therein. Recently, Pujo-Menjouet et al. [HUH] proved the existence 
of a Hopf bifurcation for the hematopoiesis model proposed in [3] , described by a nonlinear 
differential equation with discrete delay. However, their results cannot be directly applied to 
(E|) because of the nature of the delay. 

Delay differential equations with distributed delay have been studied by many authors. 
We can cite, for example, the works in [TOl ED EH EH EU • However, these studies mainly 
focused on stability conditions. In 2003, Liao et al. 15j showed the existence of a Hopf 
bifurcation for a Van der Pol equation with distributed delay and studied the stability of limit 
cycles, applying the normal form theory and the center manifold theorem. The characteristic 
equation in |15| is an exponential polynomial, similar to the one obtained in [HIH] except that 
the degree is higher, which makes the study easier than with equation (El). In [2], Adimy et 
al. obtained the existence of a Hopf bifurcation for a nonlinear differential equation with a 
delay distributed according to a density, generalizing equation (JTJ) . However, these authors 
did not study the limit cycles of their model. 

We consider a pluripotent hematopoietic stem cells population density x(t), satisfying 
the nonlinear delay differential equation (E]). The constant S accounts for natural mortality 
and cellular differentiation. The nonlinear term f3(x(t)) represents a rate of introduction 
in the proliferating phase. The last term appears to describe the amount of cells due to 
cell division: cells are assumed to divide uniformly on an interval (0,t), with r > 0, and 
dividing cells are in fact cells introduced in the proliferating phase one generation earlier. 
The assumption on the cell division comes from the fact that, even though only a little is 
known about phenomena involved in hematopoiesis, there are strong evidences (see Bradford 
et al. ES]) indicating that cells do not divide at the same age. The factor 2 describes the 
division of each mother cell in two daughter cells. 

The rate of reintroduction in the proliferating compartment j3 is taken to be a monotone 
and decreasing Hill function, given by 



in 



0" + x 



.71 



for x > 0. 



(2) 
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The coefficient (3q > is the maximum rate of reintroduction, 9 > is the resting phase 
population density for which the rate of re-entry (5 attains its maximum rate of change with 
respect to the resting phase population, and n > describes the sensitivity of (3 with changes 
in the population. This function was firstly used in hematopoiesis models by Mackey [I] in 
1978. 

In [H[S] Mackey gave values of the above parameters for a normal human body production. 
These values are 

5 = 0.05 d" 1 and = 1.77 d -1 . (3) 

Usually, n is close to 1, but Mackey [HE] reports values of n around 3 in abnormal situations. 

The value of 9 is usually 9 — 1.62 x 10 8 cells/kg. However, since we shall study the 
qualitative behavior of the pluripotent stem cells population, the value of 9 is not really 
important and setting the scale change 

*(*) - ^ 

in HI), with the function (3 given by ©, we obtain 

^L(t) = -6x(t)-(3 f(x(t)) + ^ J f(x(t+s))ds (4) 

with 

f( x ) = x>{). (5) 

1 w 1 + x n ~ w 

However, we mention that this special form of / will not be used in the following, except in 
computations in Section [4] We only assume that / is differentiable with /(0) = and, for 
x > 0, f(x)/x is decreasing and satisfies 

hm M = . 

x^ + oo X 

An illustration of such a function is presented in Fig. [TJ 

Notice that equation ((4|) has at most two equilibria: the trivial equilibrium x = and a 
nontrivial positive equilibrium x = x* . The trivial equilibrium always exists and corresponds 
to the extinction of the population. 

The nontrivial equilibrium exists if and only if 

< S < /3 /'(0) (6) 

and is then uniquely defined by 

5x* = M{x*). (7) 

This can be easily shown by using the fact that the function f(x)/x is decreasing on [0, +oo). 

Our aim in this work is to show that the unique nontrivial equilibrium of equation (JT|) 
undergoes, in a particular case, a unique Hopf bifurcation and to show the stability of the 
limit cycles following the approach in [TTl [18] . 

The paper is organized as follows. In Section [21 we establish some local stability results 
for the unique nontrivial equilibrium of ([1]) and prove that it undergoes a Hopf bifurcation. 
We study the stability of the limit cycles obtained at the bifurcation in Section [3] Our results 
are illustrated numerically in Section [H We conclude with a discussion in Section 



3 



M. Adimy, F. Crauste, A. Halanay, M. Neam^u, D. Opri§ 



Stability of Limit Cycles 




Figure 1: Graph of the function / given by ([5|) for different values of n. 



2 Local Stability and Hopf Bifurcation Analysis 

Part of the study presented in this section has been previously performed by Adimy et al. in 
[Tj. However, for the reader convenience and to preserve the coherence of the present work, 
we detail the asymptotic behavior study of the nontrivial equilibrium x = x* of equation 
dH) , defined by ([7]) . We first concentrate on the local asymptotic stability of this equilibrium. 
Then, we will show that it undergoes a Hopf bifurcation for some critical value of the time 
delay. 

We assume that (O holds in order to ensure the existence of x* , that is 

< S < o f(O). 

The linearization of equation Q around x* leads to the characteristic equation 



A(A,r) := A + 5 + o 0x - 



20001 '"' 



e Xs ds = 0, 



(8) 



where we have set 



Pi ■■= /V)- 



We recall that the nontrivial equilibrium x* is locally asymptotically stable if and only if all 
eigenvalues of (O have negative real parts. 

The function / in ^ is not necessarily monotone so 0\ may be either positive or negative. 
We first study the case 0\ > 0. 

Consider A(A,r) as a function of real A. Then A is differentiable with respect to A and 



OA 



A A (A,r) :=— (A, r) = 1 



20001 f° x 



se As ds. 



(9) 
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Since / se Xs ds < 0, we deduce that A\(\,t) > 0. Moreover, one can easily check that 



lim A(A, r) = — oo 



and 



lim A(A, r) = +00. 



A— *■— 00 A- 

Consequently, A(A,r) has a unique real eigenvalue, namely Ao- From ([5)1, we obtain 

A(0,r) =5-00/3!. 

Writing 



Pi = x* 



fix*) 



we deduce 



A(0,r) = -p x* 



> 0. 



Hence, Ao is strictly negative. 

Let A = v + iuj be an eigenvalue of (JHJ) and assume that v > Xq. Considering the real part 
of ((El), we obtain 



v + S + PoPi 



2PoPi 



e vs cos(ujs)ds = 0. 



Therefore, 



v - A 



[e vs cos(ws) - e AoS ] ds < 0, 



which gives a contradiction. We deduce that every eigenvalue of ([5]) has negative real part. 
It follows that the nontrivial positive equilibrium x* is locally asymptotically stable. This 
result is summed up in the following proposition. 

Proposition 2.1. Assume that Pi > 0. Then the nontrivial equilibrium x = x* of |7p is 
locally asymptotically stable for all r > 0. 

If the function / is given by (0 then © is equivalent to < S < p and x* — (Po/S—l) 1 ^ 1 . 
The condition /3i > reduces to 

Po-S 



Po 



< 1. 



We assume now that P% < 0. We are going to show that the equilibrium x* undergoes a 
Hopf bifurcation. To that aim, we look for the existence of purely imaginary roots of ([HJ). 

We first check that x* is locally asymptotically stable when r = 0. In this case, the 
characteristic equation {SI reduces to 



A + 5 - PoPi = 0, 



so 

A = -S + PoPi < 0. 

We have the following lemma. 

Lemma 2.1. Assume that Pi < 0. Then the nontrivial equilibrium x = x* of is locally 
asymptotically stable when r = 0. 
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Let A = iuj, u> G R, be a purely imaginary eigenvalue of (jSJ). One can check that 

A(-iw,r) =0, 

so we only look for positive o>. Moreover, w/0 since 

A(0,r) =<5-A)/3i >0. 

Thus, let us assume that lj > and t > satisfy A(ia>, r) = 0. Separating real and imaginary 
parts of A(kj,t), we obtain 



d + /3 Pi sinlUT = 0, 

LOT 

H (1 — cos(wr)) = 0. 

We set 

h(x) = — — , x > 0. 



(10) 



Then system (flTfl) can be written 



h(uT) 



2/3 /3i 

cos(wt) — 1 1 
(Wp = 2/3 /?ir" 

Since /3i < and <5 > 0, then 

g + l 

2/3o/?i < 2" 

Let 

xq = min{x > 0; x — tan(a;) and h(x) > 0} ~ 7.725, 

and assume that 



(11) 



ft(a;o) < w (12) 



One can check that h(xo) ~ 0.1284. Then (fT2| is equivalent to 

5 



Pi < 



'po(i-2h(x )y 



On the interval [0, 7r], the function /i is strictly decreasing and nonnegative, with < h(x) < 1. 
Moreover, for x > 7T, < h(xo)- Consequently, the equation 

Mx) = 

has a unique solution, denoted x c , which belongs to the interval (0, tt). We then set 



2/3o/3i(cos(a; c ) - 1) 

and 



U c — — 
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Therefore, (lo Ci t c ) is the unique solution of and ±iuj c are purely imaginary eigenvalues 
of © for r = T C . 

In order to show that x* undergoes a Hopf bifurcation for r = r c , we have to prove that 
±iuj c are simple eigenvalues of A(-,t c ) and satisfy the transversality condition 



dRe(A) 



dr 



> 0. 

T = T C 



We first check that ±ioj c are simple eigenvalues of (jHJ ■ Using ([9]) , one can see that ico c is 
simple if 

t} /a c \\ i , uj c T c sm(u c T c ) + cos(cj e T c ) - 1 

Re(A A («w c ,T c )) = H j ^° ( 13 ) 

or 

Im(A A (iw c ,T c )) = ^ / 0. (14) 

We are going to show that, in fact, these two conditions are satisfied. 

Lemma 2.2. Assume that f3\ < and dig)) holds. Let (uj c ,t c ) be the unique solution of hll}) . 
with lo c t c E (0,7r). TTien 

i?e(A A (iw c ,T c )) > and Im (A x (iuj c , r c )) > 0. (15) 

/n particularly, ±iuj c are simple eigenvalues of ^ /or r = r c . 

Proof. First, one can check, using (fTT]) and (fT5|) . that 

Re (A a (zu; c ,t c )) = 2 + (5 + /3 /3i)r c . 

Since 

cos(ir) — 1 h 2 (x) 

H = -t i t \ i for * > 0, 

ar 1 + cos(x) 

then, from (fTT|) it follows that 

1 + cos(x) = -(5 + PoPi)t c - 

Consequently, 

Re (A A (iw c , t c )) = 1 - cos(cj c t c ) > 0. 

Secondly, since xcos(x) < sin(x) for x E (0,7r) and x c — ui c t c E (0,7r), then from (jT^J) we 
obtain 

Im(A A (^ c ,T c )) > 0. 

This concludes the proof. □ 

Consider now a branch of eigenvalues A(r) = v(t) + iuj(r) of ([5]) such that v{t c ) = and 
lo(t c ) — uj c . Separating real and imaginary parts in ^ we obtain 

v{t) +S + /3 f3 1 - J° e"( T > cos(u(T)s)ds = 0, 

u {r)-^Ml f e v ^ s sm{oj( T )s)ds = 0. 
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Then, by differentiating each of the above equalities with respect to r, we get, for r = t c , 
Be(Ax(iu C) T c ))u'(T c ) 

t /a c \\ a \ , 2 PoP\ ( , , sin(a; c ) 
= Im(A A (iw c ,T c ))o/(T c ) H cos(x c ) 



(16) 



and 



2/3 /3i (I- cos(x c ) 



Re(A A (iw c ,r c ))w'(T c ) 
= -Im(A A (wj c ,T c )) v'(t c ) + 



Using (TTK|) . (Unj) and (TT7|) . we can see that f'(r c ) satisfies 



Im(A A (icj c ,T c )) 2 + Re (A A (?o; c , t c )) 2 i/(t c ) 



— sin(x c ) 



(17) 



1 — cos(a; c ) 



cos(x c ) 



sin(x c )J Im(A A (iw c ,T c )) 
»Re(A A (io; c ,T c )) 



Using the definitions in (fT3|) and (I14|) . simple computations give 
1 — cos(x c ) 



+ cos(a; c ) - 



sin(x c ) ) Im(A A (iw c ,T c )) 
sin(a; c ) 



Re (A\(iu c ,T c )) = cos(x c ) - 



sin(x c ) 



Hence, 



It follows that 



a /. n,2 // x 2/3 /3i x c cos(x c ) - sin(a; c ) 
A A («w c , t c )\ v (r c ) = > 0. 



v'{t c ) > 0. 



(18) 

To conclude, when r = r c , the characteristic equation A(A, r) has a unique pair of purely 
imaginary simple eigenvalues satisfying (dRe(A)/dr)(r = r c ) > 0. Consequently, a Hopf 
bifurcation occurs at x* when r = r c . Moreover, applying Rouche's Theorem with Lemma 
12.11 we easily check that every eigenvalue of A(A,r), with r < r c , has negative real part. It 
follows that x* is locally asymptotically stable for < r < r c . These results are summed up 
in the following theorem. 



Theorem 2.1. Assume that flx < and hlty) holds. Then there exists a unique value r c > 
of the time delay such that the equilibrium x = x* is locally asymptotically stable when 
t G [0, t c ) and becomes unstable when r = r c throughout a Hopf bifurcation. In particularly, 
periodic solutions appear for equation |7p when r = t c . 

As an example, one can check that when / is given by ([5]) the assumptions in Theorem 
12.11 are equivalent to 



n > 



2(1 - h(x )) fa 



2.35- 



l-2h(x ) O -S ' Po-8' 

In particularly, these conditions are satisfied when (3q and 5 are given by ([3|) and n > 2.42. 

The existence of a Hopf bifurcation in Theorem 1 2 . 1 1 leads to the existence of a limit cycle 
when the bifurcation occurs. In the next section, we focus on the stability of this limit cycle. 
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3 Stability of Limit Cycles 

We study now the stability of the limit cycle yielded by Theorem l2.ll We follow the approach 
used in [TTJITH] - This involves the description of a center manifold and subsequently the study 
of the normal form given by the restriction of the flow to this center manifold. The stability 
of the limit cycle will be decided by the sign of the first Lyapunov coefficient l\ (0) . 

For general properties concerning delay equations and the theory of central manifolds for 
these equations, see |19j . For the existence and various properties of center manifolds we 
refer to [HI El H2 H3 US US] • Also, in [53] and [H], a rigorous treatment of the operators 
involved in this approach is to be found. A similar problem is considered in [15 . 

Define, for t > 0, 

y(t) = x(t) - x*, fi = T-T c 

with x* the nontrivial equilibrium of Q that bifurcates into a limit cycle for the critical 
value t — t c (see Theorem 12. 1|) . The equilibrium x* is defined by ([6|) and ((7|). Equation (j4|) 
turns into 

y\t) = -S(y(t) + x*) - (3 f(y(t) + x*) + ( f (y(t + s) + x*) ds. (19) 

Thanks to this formulation, we now concentrate on the trivial equilibrium y = of (|19p 
which bifurcates when fi = 0. 

For an interval I C K, denote C(J,K) = {/:/—»■ K, / continuous } where K = K or C. 
When I = [—fj, — t c , 0], we set 

G p :=C([-m-t c ,0],K). 

Considering, for t > 0, the function y t : [—fj, — r c , 0] — > K defined by yt(s) = y(t + s), we 
can reformulate equation (|19p as the following abstract functionnal differential equation 

j t y{t) = GM, t > 0, (20) 

where, for tp £ C^, 

G^y) = -.5 [^(0) + x*} - fof (cp(0) + x*) + [ f (<p(s) + x*) ds. 

Assume that / is C 4 on [0, +oo) (remark that G^ is then C 4 (C At ,K)). 

Consider the linearized equation of ([20|. corresponding to the Frechet derivative D v G M (0) := 
Lfj,, given by 

j t z(t) = L^Zt, t > 0. (21) 

In fact, is given explicitly by 

= -cMO) + c 2 (/x) / ¥>(0)d0, 9? e G M , (22) 

J —fX—T c 



where 



Setting 



Cl :=<5 + /? /3i, c 2 (/x):=^^, ft := f'{x*). (23) 
Ffi '■= G M — 
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equation P0|) becomes 

-y(t) = L ll y t + F ll (yt), t>0, (24) 

with F M (0) = and D^F^O) = 0. 

In order to develop a normal form associated to equation (|2ip , we write this latter as an 
abstract ordinary differential equation. 

First, we know from [24] that the linear equation (f2T|) gives a Co-semigroup (T(t)) t >o on 
Cp, with generator A^ defined by 

( V(A^) = {pgCHl-M— r c ,0],R); v '(0) = L^} , 



To write equation (|24|) as an ODE we need to extend the problem (|21() to the Banach space 
C,, := 8 (X ), where 

(X ) = {X c; eel and (X c)(6») = X ((9)c} 

and Xq denotes the function defined on [— ^ — r c , 0] by 



*o(0) = 



0, if - (i - t c < e < o, 

1, if = 0. 



Adimy proved in [37] that this extension determines a Hille-Yosida operator. This result is 
recalled in the next lemma. 

Lemma 3.1. The continuous extension A^ of the operator A^ defined on C M by 

(25) 

A^ = ip' + X (L li <p-<p'(0)), <p£V(A fi ), 

is a Hille-Yosida operator on C M ; that is: there exists luq £ K suc/i t/iai ((x>o,+oo) C p(^) 
and 

sup | (A - u ) n \\(XI - \)- n \\,n e N, A > w } < oo. 

It follows that if y is a solution of (j24|) on [0, T], T > 0, with an initial condition (p £ 
on the interval [— /x — r c , 0], then the function t £ [0, T] i— > yt G C M satisfies 



d 

Jt"" ' " ^ (26) 



-y t = A^yt+XoF^yt), t £ [0,T], 



Conversely, if there exists a function t £ [0, T] \— > G C M such that 

— (t) = ^u(t)+X F p (u(t)), ie[0,T], 



u(0) = 
then = j/t, t S [0, T], where 



y(*) 



u(t)(0), if*€[0,n 
¥»(*), if * e [-At-r c ,0], 
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and y is a solution of (|24p . This handles in particularly the problems arising from the fact 
that Ajj, does not preserve the space of ^-functions. 

Thanks to results by Arendt [25J and Da Prato and Sinestrari [29 , the ODE (J2TJ) is 

well-posed for initial conditions in T>(A^) = C\ t . 

Now we can reformulate the problem (|24[) as the abstract ODE ([2"o]) . 

Another important step towards the description of a center manifold is the definition of 
a bilinear form related to the equation (j2"TT) . 

From now on, we set IK = C. For cp G C M and ip G C* := C([0, /j, + t c ], C), define according 
to [23] or [24], 

r 



<tM = ^(omo)- 



/i-T c 



where d?7(s) = c 2 (fi)ds ~ c\X§(s). Thus ([25]) becomes 



(28) 



(29) 



We build a natural extension of this bilinear form to the space C* x C M where 

c; = c; © (x *) 

and 

(X *) = {X*c; c G C and (X*c)(9) = X*(9)c} 
with Xq the function defined on [0, [i + t c ] by 

f 0, if < 6< fi + t c , 
ow \ 1, if 6 = 0. 
We obtain, for tp € C*, ip £ and a, c 6 C, 

(V> + Xq a, </? + A c) = (V>, <p) + ac. 
With respect to this bilinear form, we define the adjoint of the operator A M , denoted A* 
and its domain V(A*). It satisfies, for tp G £>(A M ) = ^([-/x - t c , 0], C) and ip G 2>(A*), 

From ([25]) and ([29]). we obtain 



(-!/>, A^i/j) = ip(0)L^(p - c 2 (m) 



Using an integration by parts and (|22|) , we deduce 



(tM/^) = ^(0) 



-ci<^(0) + c 2 (/x) 



¥>(0)d0 



/i-T c 



-C2(M) 



/i-T c 



^(0)y(a) -V(-sMO)- 



^V(£-sM£R) da, 



-ci^(0) + c 2 (/x) 



At— r c 



p(0) 



+c 2 (^) 
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where 




&([0,H + T e ], 



C2(M) 



^(s)ds-ci^(0) + V'(0) 



We consider now the purely imaginary eigenvalues of jSJ denoted ±iuj c , with oj c > 0, 
which exist when r = r c , that means when the bifurcation occurs (see Section [5] and, in 
particularly, Theorem 12. ip . From the definition of the characteristic equation in ((HI) and the 
notations introduced in (|23|) . we have 



It follows that 



A(«w c , r c ) = iw c + ci - c 2 (0) 

c 2 (o) [ - L M- 



s ds = 0. 



ci = ioj c . (30) 
G C 1 ([— r c , 0], C) is an eigenvector of A 



Then, with definition (|25l) . the function q(s) = e 4 ' 
associated with iu c . 

Hence, q*(s) — de lulcS £ C 1 ([0, r c ], C), d 7^ 0, is an eigenvector for Aq associated with 
—ilo c . Moreover, we can choose d G C so that the norming condition (q* , q) = 1 is satisfied. 
It follows that 



d = 



1 + 02(0) 



One can check that in fact d — (A\(iuj c ,T c ))~ . Since ilu c is a simple root of A(-,t c ), then d 
is well-defined. 

From (l29l) and (1301) we infer also that 



(<?*,?> =0. 



(31) 



We are interested in the center manifold corresponding to the eigenvalue A = ilu c of Aq and 
to the system (|26|) . Such a center manifold exists (see [HI [21]): it is a locally invariant, locally 
attracting manifold containing the origin and tangent at the origin to the subspace spanned 
by the eigenvectors corresponding to the eigenvalues ±itu c of Aq. In fact, to reach our aim, 
we only need information on the section of the center manifold, denoted Co, corresponding 
to n = (see [T7]). 

Let yt be a solution of 

^=A yt + XoFo(y t ). (32) 

We compute the coordinates of the section Co of the center manifold corresponding to fi = 0. 
Following the notations in IT] , we define 



z (t) = (q*,yt), 



for t > 0. 



(33) 



We will use z and z as local coordinates of Co in the directions q* and q* respectively. We 
also define, for t > and s G [— t c , 0], 

w(t, s) = y t (s) - z(t)q(s) - z(t)q(s), 
= y t (s)-2Rc[z(t)q(s)}. 



We have 



w(t,s)=W(z{t),z(t),s), t>0, sg[-t c ,0], 
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with 



z 2 _ z 2 

W{z, z, s) = w 20 (s) — + wu(s)zz + w 02 (s) — + 



(34) 



One can notice that, for real solution y, w is real so wq2 — wio- Moreover, (|31j) and (|33p 
imply that (q*,w) = 0. 

The section Co of the center manifold is locally invariant under equation (|32j) : any solution 
that starts in it will stay in it for any time t in some nontrivial interval; therefore, if y t G Cq 
we have 

^-z(t) = (q*,A y t +X F a (y t )) 



dt 



so, from 



it follows that 
d 



dt 



z{t) = iuj c z(t) + dF (W(z(t),z(t), •) + 2Re[z(t)q}) , 
= iu c z(t) + g(z(t),z(t)), 



(35) 



with 

g{z,z) = dF a (W{z,z, •) + 2Re[zg]) . 
We use the Taylor expansion of / around a;* to rewrite Fq as 



where 



Fo(tp) = c 3 



v{e) 2 



3! 



tp(6f + G(<p(d) 4 



21 



^(0) 2 



3! 



^(0) 3 + O((p(0) 4 ) 



2/?o 



dB 



c 3 := — , fa := /"(x*) 



and 



A := /"'(a:*). 



If we denote, for convenience, w(s) = W(z, z, s), we then obtain 



(36) 



g[z,z) = rfc 3 <{ y / [w(s) + ze^" 8 + ze- luJ " s ] 2 ds 



6 



[w(s) + ze w " s + ze-^ s ] 3 ds 



iuj c s i —„—ioj c s~\4: 



)ds 



(37) 



0([tu(s) +ze 

- p o d(^-[w(0) + z + z} 2 

+ y[w(0)+z + z] 3 + O([ W (0) + z + z] 4 ) 

Equation |33|) is called the normal form obtained by the restriction of the flow to the 
center manifold. Our next goal is to compute some coefficients in the Taylor series of g and 
to use them to study stability of the limit cycle by computing also the Lyapunov coefficient. 
This latter is given by some coefficients in the Taylor expansion of g(z,z) given by (|37j) . This 
means, in fact, that stability of the limit cycle is determined by the normal form obtained 
through the restriction of the flow to the center manifold. 
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Restricting the Taylor expansion in (|37[) to terms of order less or equal to three, we get 
g( z , z) = ^g 20 z 2 + gxxzz + ^g 02 z 2 + \^g 2 \z 2 z + ... (38) 



with 



320 = -d(3 2 (fo + C3(1 ^P^- i 
gn =df3 2 {c3T c - A))i 

go2 = -dfh (A ^ l j , m 

92i=d(c 3 I [ [w 20 ( S )e^ s + 2w 11 ( S )e l ^ s }ds - /3 3 (1 ~ 6 ^ 1 



-/3o/32[w 20 (0) + 2 Wll (0)]-/3 /3 3 
With the definition of C3 given in (|36[) . one can see that 

.911 = d0 o 2 - 

One can note that we do not need the coefficients gij with i + j > 2 except 321 (see [T7]). 

At this point, we still need to compute w 2 q(s) and wn(s) for s G [— r c ,0]. It follows 
directly from (J32J), ((33D and (J33J), that 

!»(*,-) = f t y t -§ t [z(t)q-z(t)Q], 



(40) 



= A w(t, •) + X F (w(t, •) + 2Re(z(t)q)) 

-2Re (g(z(t),z{t))q), 
= A w(t,-) + H(z(t),z(t),-), 

where, for s G [— r c ,0], 

H(z,z,s) = -2Re(g(z,z)q(s))+X (s)F (W(z,z,-) + 2Re(zq)). (41) 
Let s G [-r c ,0) be fixed. From (gTJ) with dHHJ> we have 
H(z,z,s) = -2Re (g(z,z)q(s)) , 

z 2 z 2 \ 

520y + gnzz + g Q2 — + ...J q{s) 

z 2 _ _ _ z 2 \_ 
#2oy +3n zz + .9o2y + •••) q( s )- ( 42 ) 



Considering the expansion 



H(z,z,s) = H 2 o{s)—+H n {s)zz + H 02 {s)— + 



and comparing the coefficients with those in (|42f yields 

H 20 (s) = -g 2 oq{s) -g 02 q(s), 

Hu(s) = -gnq(s) - g n q{s), (43) 
H 02 (s) =H 20 (s). 



14 



M. Adimy, F. Crauste, A. Halanay, M. Neam^u, D. Opri§ 



Stability of Limit Cycles 



From (O, ([25]) and |@0J), we obtain 

(^4 - 2iu c )w 2 o{s) = -H 2 o(s), 

A w n (s) = -H u (s), (44) 
(A + 2iw c )w 2(s) = -H 02 (s). 
Using and gSj), flUD becomes 



ds 
dwi 



ds 

Solving the above system, we obtain 



■(a) = 2icj c w 2 o(s) + g2og(s) + g 2<l( s ), 
-(s) = ffiig(s) +5 n ?(s). 



W20 ( s ) = _p9. e *-.» _ J^Lg-w + Ei^""', (45) 
«, 11 (a) = ^l e *-o._lii e -fa, . + _ B (4g) 



where E\ and i?2 can be determined by setting s = in (|4ip . In fact, we have 
iJ(z, z,0) = -2Re (g(z,z)) + F (W(z,z, •) + 2Re(zg)) , 



so we deduce 

#20 (0) = -.920 - ff 02 + Ih [ c 3 / e^'ds - #, ) , 





2iuj c s , 



Hn (0) = - g n + fa ( c 3 / ds-(3 



From (I25l) and (HI), we have 



c 2 (0) ( w 2Q {s)ds-{c x +2iu c )w 2Q {Q) = -H 20 (0), (47) 

J—Tc 

cb(0) / «;ii(s)ds-ci«;ii(0) - -#n(0). (48) 

Substituting (|4*5|) in (|4"T|) and using (I5U|) we eventually get 

R 02(c 3 (l-e-^)-2tu;M 
1 c 2 (0)(l - e- 2<w ^) - 2ciio; c + 4^2 ' 1 ' 

Similarly, substituting (|4*U]) in (|4*5|) we get 

£ 2 = - Mffi^l . (50) 
c 2 (0)r c - ci 

Using the definitions given in (f2Uf and (|5S|) . we can write 

/3o/3 2 
2 ^ - 0o0i ' 
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We are now able to complete the calculation of 521 in (J2SJ) using the above values of w 2 o, 
wn, Ei and E 2 . We set 

r° 

K:= \ [w 20 {s)e-^ s + 2w n (s)e^ s ]ds. 



Then, from (J45J) and l[lg]). we obtain 

(g 02 +6 gil )(l-e- 2 ^) r e ( g20 + 2g 11 )-(£; 1 + 2£; 2 )(l-e-^) . 
A = — 1 1 (51) 

where E\ and E 2 are given respectively by |49|) and (|50l) . 
From pg]l. IjlSjl and (gBJ, we have 

521 = d c 3 /3 2 if - c 3 /3 3 i 

a a ( t? lop 1 3520 + 502 • , 2 (5n ~ 511 ) A o „ 
-PoP 2 -Ei + 2£? 2 H 7, 1 H * - PoA 

where if is given by (f5Tj) . 

Based on the above analysis and calculation, we can see that each gij in (|39|) is determined 
by the parameters and delay in equation Q. Thus we can explicitly compute the following 
quantities: 

£i(0) = (520511 ~ 2|gii| 2 - ^|ffo2| 2 ^ + ^521, 

li(0) = Re(Li(0)), 

^ 2 = -ReTyR)' (52) 
b 2 = 2ii(0), 

Im(ii(0))+/i 2 ImA'(r c ) 



2 



One knows from [T7] that the following properties hold: if /i 2 > (< 0) then the Hopf 
bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for r > t c 
(t < t c ); solutions are orbitally stable (unstable) if b 2 < (> 0); and the period of bifurcating 
periodic solution increases (decreases) if T 2 > (< 0). 

The coefficient A'(r c ) in ([52")) is given by (fH>)) and (jTT|> . In particularly, we have proved in 
Section property (fT5|) . that 

Re(A'(r c )) > 0. 
In summary, this leads to the following result: 

Theorem 3.1. If the Lyapunov coefficient h(0), defined in i52\) . is negative (resp. positive) 
then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions 
exist for t > t c (t < t c ), and solutions are orbitally stable (unstable); The coefficient T 2 
determines the period of the bifurcating periodic solutions: the period increases (decreases) if 
T 2 > (< 0). 
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4 Numerical results and simulations 

We numerically compute, in this section, the formulas obtained above to determine the 
behavior of the periodic solutions of equation Q. 

We choose / as in ([5]) . In order to satisfy the assumptions of Theorem 12. 1[ we have to 
choose 5, f3o and n such that (3\ < and (fT^j) holds true; that is 

2(l-h(x )) A) 
l-2h{x ) [3 -5- 

We take S and (3q as given in Then the above conditions are in fact satisfied for n > 2.42. 

Using Maple 9, we are able to compute the coefficients in (|52p . listed in the following 
table for n G {3,4,5,6,7,8}: 



n 


X* 






h(0) 


T 2 


3 


3.2523 


18.1270 


0.1380 


-0.0026 


0.1540 


4 


2.4218 


10.9688 


0.2091 


-0.0375 


0.8083 


5 


2.0291 


7.8748 


0.2785 


-0.1504 


2.2097 


6 


1.8034 


6.1447 


0.3472 


-0.3975 


4.5441 


7 


1.6577 


5.0385 


0.4155 


-0.8399 


7.9335 


8 


1.5562 


4.2702 


0.4836 


-1.5412 


12.4589 



Even though we do not give values of h(0) and Ti for all n > 2.42, we can notice that 
observations indicate that /i(0) is strictly negative and T<i strictly positive for n > 2.42. 
Hence the unique Hopf bifurcation of equation ([4]) seems to be supercritical and solutions 
orbitally stable, with increasing periods. 

Using the Matlab solver ddc23 [30 , we can compute the solutions of equation fl]) for the 
above-mentioned values of the parameter n and for any positive initial condition. Solutions 
of are shown versus time t and in the phase plane in Fig [5] to [7] 




R °° 100 200 ™ 400 500 " °° 2 4 -» 

Figure 2: Solutions of (|U) are displayed when n — 3. Periods of the oscillations are close to 
46 days. 



5 Discussion 

Many hematological diseases involve oscillations about a steady-state during the chronic pe- 
riod. These oscillations give rise to instability in the hematopoietic stem cell count. Chronic 
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Figure 4: Solutions of are displayed when n = 5. Periods of the oscillations are close to 
23 days. 




Figure 5: Solutions of ^ are displayed when n = 6. Periods of the oscillations are close to 
18 days. 
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Figure 6: Solutions of (jl]) are displayed when n = 7. Periods of the oscillations are close to 
15 days. 




Figure 7: Solutions of (j4]) are displayed when n = 8. Periods of the oscillations are close to 
13 days. 
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myelogenous leukemia (see Fortin and Mackey [3T]) is one of the most common types of 
hematological disease characterized by the existence of periodic oscillations (oscillations of 
leukocytes with periods from 30 to 100 days). Experimental observations have led to the 
conclusion that this dynamic instability is located in the hematopoietic stem cells compart- 
ment. 

We have studied, in this paper, a mathematical model of pluripotent hematopoietic stem 
cells dynamics in which the length of the proliferating phase is uniformly distributed on an 
interval. We have shown that instability can occur in this model via a Hopf bifurcation, 
leading to periodic solutions usually orbitally stable with increasing periods. This has been 
obtained throughout the description of a center manifold and the subsequently study of the 
normal form. 

Periods of the oscillations obtained in numerical simulations, in Section 31 may be in the 
order of 30 to 50 days (at the bifurcation) when the parameter n is not too large, correspond- 
ing to what can be observed with chronic myelogenous leukemia. It has already been noticed 
by Pujo-Menjouet and Mackey [5] that this parameter n, which describes the sensitivity of 
the rate of reintroduction j3, plays a crucial role in the appearance of periodic solutions when 
the delay is constant. The sensitivity n describes the way the rate of introduction in the 
proliferating phase reacts to changes in the resting phase population produced by external 
stimuli: a release of erythropoietin, for example, or the action of some growth factors. Since 
periodic hematological diseases are supposed to be due to hormonal control destabilization 
(see |31j). then n seems to be appropriate to identify causes leading to periodic solutions. 
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